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Abstract 

The  computation  of  the  analytic  center  of  the  solution  set  can  be  important  in  linear 
programming  applications  where  it  is  desirable  to  obtain  a  solution  that  is  not  near  the 
relative  boundary  of  the  solution  set.  In  this  work  we  discuss  the  effective  computation 
of  the  analytic  center  solution  by  the  use  of  primal-dual  interior-point  methods.  A 
primal-dual  interior-point  algorithm  designed  for  effectively  computing  the  analytic- 
center  solution  is  proposed  and  theory  and  numerical  results  are  presented. 


1  Introduction 

The  computation  of  the  analytic-center  of  the  solution  set  is  valuable  for  some  linear  pro¬ 
gramming  applications  where  it  is  desirable  to  obtain  a  solution  that  is  not  near  the  relative 
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boundary  of  the  solution  set  (see  Charnes,  Cooper,  and  Thrall  [3]  for  such  an  application).  In 
this  work  we  discuss  the  use  of  primal-dual  interior-point  methods  for  effectively  computing 
the  analytic-center  solution. 

Because  in  this  work  we  are  interested  in  using  primal-dual  interior-point  methods  for 
finding  a  specific  solution  (the  analytic-center  solution),  the  study  of  the  convergence  of  the 
iteration  sequence  generated  by  primal-dual  interior-point  methods  needs  to  be  addressed. 
It  is  relevant  to  mention  that  this  is  not  the  case  for  many  real-world  problems  where  any 
feasible  point  with  the  value  of  the  duality  gap  close  to  zero  is  a  satisfactory  approximate 
solution. 

The  basic  primal-dual  interior-point  method  for  linear  programming  was  proposed  by 
Kojima,  Mizuno,  and  Yoshise  [11].  Most  primal-dual  interior-point  methods  fall  into  the 
general  framework  of  the  Kojima-Mizuno- Yoshise  algorithm.  The  first  work  devoted  to  the 
study  of  the  convergence  of  the  iteration  sequence  generated  by  the  Kojima-Mizuno- Yoshise 
algorithm  is  Tapia,  Zhang,  and  Ye  [18].  Later  on,  Zhang  and  Tapia  [24]  derived  conditions 
under  which  this  iteration  sequence  converged  to  the  analytic-center  solution,  assuming  that 
the  sequence  converged.  However,  the  conditions  given  in  [24]  are  not  compatible  with  the 
Tapia-Zhang- Ye  conditions  for  the  convergence  of  the  iteration  sequence.  The  former  requires 
the  centering  parameters  to  be  bounded  away  from  zero,  while  the  latter  requires  that  these 
same  parameters  converge  to  zero. 

An  interesting  variant  of  the  Kojima-Mizuno- Yoshise  algorithm  was  proposed  by  Mizuno, 
Todd,  and  Ye  [16].  The  Mizuno- Todd- Ye  predictor-corrector  algorithm  is  a  path  following 
algorithm  which  keeps  the  iterates  in  a  small  neighborhood  of  the  central  path.  Zhang  and 
El-Bakry  [22]  showed  that  a  modified  version  of  the  Mizuno- Todd- Ye  predictor-corrector 
algorithm  had  the  property  that  the  iteration  sequence  that  it  generated  converged  to  the 
analytic-center.  Soon  after  the  result  of  Zhang  and  El-Bakry,  Gonzaga  and  Tapia  [8]  proved 
that  the  Mizuno- Todd- Ye  predictor-corrector  algorithm  itself,  without  any  modification,  gen¬ 
erated  an  iteration  sequence  which  converged  to  the  analytic-center.  These  are  the  first 
results  proving  convergence  of  the  iteration  sequence  to  the  analytic  center  solution  for  any 
primal-dual  interior-point  method. 

The  Mizuno- Todd- Ye  predictor-corrector  algorithm  has  two  characteristics  that  makes  it 
unappealing  for  most  practical  purposes.  The  first  one  is  that  at  each  iteration  the  algorithm 
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requires  the  solution  of  two  linear  systems  of  equations  where  the  matrices  are  of  the  same 
form  as  those  for  the  Kojima-Mizuno-Yoshise  algorithm.  Hence,  the  cost  per  iteration  is  twice 
that  of  the  Kojima-Mizuno-Yoshise  algorithm.  The  second  deficiency  is  the  requirement  that 
the  starting  feasible  point  must  be  inside  a  specified  /^-neighborhood  of  the  central  path. 

Recently,  Bonnans  and  Gonzaga  [1]  gave  sufficient  conditions  for  convergence  and  con¬ 
vergence  to  the  analytic-center  of  the  iteration  sequence  generated  by  primal-dual  interior- 
point  algorithms  from  a  general  class  of  algorithms  containing  the  Kojima-Mizuno-Yoshise 
and  Mizuno-Todd-Ye  algorithms.  Their  result  is  an  extension  of  the  work  of  Gonzaga  and 
Tapia  [8].  We  have  found  that  theoretical  results  are  often  limited  in  practice.  Theorem  3.2 
in  Bonnans  and  Gonzaga  [1]  guarantees  the  theoretical  convergence  of  the  iteration  sequence 
to  the  analytic-center  assuming,  among  other  conditions,  that  the  iterates  stay  in  a  small 
neighborhood  of  the  central  path.  However,  strong  adherence  of  the  iterates  to  the  central 
path  results  in  slow  performance  of  the  algorithm,  especially  when  this  adherence  is  also 
maintained  far  from  the  solution  set.  In  addition,  we  will  show  that  although  the  iteration 
sequence  generated  by  the  Kojima-Mizuno-Yoshise  algorithm  may  theoretically  converge  to 
the  analytic-center,  in  practice  this  convergence  is  not,  for  many  cases,  easy  to  attain.  This 
negative  behavior  is  related  to  the  proximity  of  the  iterates  to  the  solution  set.  If  the  se¬ 
quence  generated  by  the  algorithm  approaches  the  solution  set  very  fast  and  is  far  from  the 
analytic-center  the  stopping  criteria  are  met  before  the  approximate  solution  has  a  chance 
of  getting  close  to  the  analytic-center.  Also,  the  singularity  of  the  Jacobian  matrix  at  the 
solution  may  cause  serious  numerical  problems  that  affect  the  computation  of  the  analytic 
center  when  the  iterates  are  very  close  to  the  solution  set. 

Hence,  a  balance  between  proximity  to  the  central  path  and  proximity  to  the  solution 
set  is  crucial  for  an  effective  computation  of  the  analytic-center  of  the  solution  set  using  a 
primal-dual  interior-point  algorithm. 

In  this  work  a  modification  of  the  Kojima-Mizuno-Yoshise  primal-dual  algorithm  is  pro¬ 
posed.  This  modified  algorithm  combines  the  objectives  of  approaching  the  central  path 
and  the  solution  set  to  effectively  compute  the  analytic-center  solution.  We  prove  that  the 
proposed  algorithm  generates  an  iteration  sequence  with  the  property  that  a  subsequence 
converges  to  the  analytic-center  of  the  solution  set;  however  in  practice  we  always  have  that 
the  entire  sequence  converges  to  the  analytic  center.  Moreover,  we  also  prove  that  with  a 
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small  modification  the  algorithm  achieves  polynomial  complexity.  Finally,  we  demonstrate 
that  the  numerical  behavior  of  the  algorithm  is  quite  good. 

This  paper  is  organized  as  follows.  The  next  section  contains  some  preliminary  back¬ 
ground  and  notation.  In  Section  3  our  proposed  primal-dual  interior-point  algorithm  de¬ 
signed  for  effectively  computing  the  analytic-center  is  described.  The  convergence  properties 
for  the  algorithm  are  studied  in  Section  4.  In  Section  5  we  prove  that  the  algorithm  has  poly¬ 
nomial  complexity  if  the  algorithm  is  modified  in  an  appropriate  way.  In  Section  6  we  study 
the  numerical  behavior  of  the  proposed  algorithm.  Also,  we  compare  the  numerical  behav¬ 
ior  of  the  proposed  algorithm  with  the  numerical  behavior  of  the  Kojima-Mizuno-Yoshise 
and  the  Mizuno-Todd-Ye  algorithms.  Finally,  some  conclusions  and  remarks  are  given  in 
Section  7. 

2  Preliminaries  and  Notation 

We  consider  the  linear  programming  problem  in  the  standard  form 

minimize  cTx 
subject  to  At  =  b,  x  >  0, 

where  c,  x  £  Rn,  b  £  Rrn,  A  £  RmXn(m  <  n). 

The  Karush-Kuhn-Tucker  optimality  conditions  for  (1)  are 

*  Ax  —  b  ^ 

F{x,y,z)=  AT y  +  z  —  c  =0,  (x,z)>0, 

v  XZe  j 

where  X  =  diag(x),  Z  —  diag(z ),  and  e  =  (1, 1,  ...1, 1)T  £  R". 

The  feasibility  set  of  problem  (2)  is 

T  =  {(x,y,z)  :  Ax  =  b,ATy  +  z  =  c,  (x ,z)  >  0}. 

A  feasible  point  (x,y,  z)  £  T  is  said  to  be  strictly  feasible  if  x  and  2  are  strictly  positive. 
The  set  of  all  the  strictly  feasible  points  is  denoted  by  JF+. 


(1) 


(2) 
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We  denote  the  solution  set  of  problem  (2)  by 

S  =  {( x,y,z )  :  F(x,y,z)  -  0 ,(x,z)  >  0}. 

A  point  ( x,y,z )  in  <S  is  said  to  be  a  strict  complementarity  solution  if  x,  +  Z{  >  0  for  all 
i  —  1, . . . ,  n. 

The  following  standard  conditions  are  assumed: 

(Al)  {.t  6  Rn  :  Ax  =  b,  x  >  0}  ^  0,  and 
(A2)  Uy,z)  e  Rm+n  :  ATy  +  z  =  c,z>  o}  /  0. 

(AS)  A  has  full  rank  m. 

We  are  particularly  concerned  with  the  case  when  S  is  not  a  singleton  set. 

Under  the  previous  assumptions,  the  solution  set  S  has  the  following  interesting  structure: 

(i)  The  set  S  ^  0  is  bounded. 

(ii)  All  points  in  the  relative  interior  of  S  are  strict  complementarity  solutions  and  all 
points  on  the  relative  boundary  of  S  are  not. 

(iii)  The  zero-nonzero  pattern  of  points  in  the  relative  interior  of  S  is  invariant. 

See  El-Bakry,  Tapia,  and  Zhang  [6]  for  proofs. 

Therefore,  for  any  (x*,y*,z*)  £  ri(S),  where  ri(S )  denotes  the  relative  interior  of  S ,  the 
following  index  sets 

/+  =  {i  :  x*  >  0, 1  <  i  <  n}  and  If  =  {i  :  z*  >  0, 1  <  i  <  n} 
are  independent  of  the  choice  of  (x*,y*,z*).  Moreover,  by  strict  complementarity 

iJU'M1.2 . »)  “d  4+fU+  =  »- 

Because  of  this  structure  of  the  solution  set  S ,  the  solutions  in  ri(S)  may  be  characterized 
as  having  a  maximal  number  of  nonzero  components.  Among  them,  there  is  a  solution  that 
may  be  thought  of  as  the  center-most  solution  in  the  sense  that  it  maximizes  the  product  of 
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the  positive  components.  This  solution  is  called  the  analytic-center  of  the  solution  set  and 
was  studied  by  McLinden  [13]  in  a  general  setting  and  later  independently  by  Sonnevend 
[17]  in  the  context  of  linear  programming.  Formally,  the  analytic-center  of  the  solution  set 
S  is  defined  as: 

( x*,y*,z *)  =  arg max  {ip(x,z)  :  ( x,y,z )  G  S }  (3) 

where 

Hx,z)  =  n  x«  n  2«- 

iat  ieit 

Equivalently,  in  (3)  one  can  replace  ij)(x,z)  by  its  logarithm,  i.e., 

In  ij)(x,z)  =  E  lnx*  +  E  ln^- 

tel}  iei } 

The  central  path  associated  with  problem  (2)  parametrized  by  the  parameter  y  is  defined  as 
the  collection  of  points  Vc  =  {Vc(y)]  y  >  0}  where 

Vc{n)  —  (x(y),y(y),z(y))  is  a  strictly  feasible  point  satisfying  X(y)Z(y)e  =  ye. 

This  is  equivalent  to  saying  that  a  strictly  feasible  point  ( x,y,z )  is  on  the  central  path  (for 
some  y  >  0)  if  and  only  if  it  satisfies  xiZi  =  x2z2  =  ...  =  xnzn. 

These  notions  of  analytic-center  and  central  path  are  well-defined  under  assumptions 
(Al),  and  (A2).  For  more  details  see  McLinden  [13],  Megiddo  [15],  and  Gonzaga  [7]. 

A  /^-neighborhood  of  the  central  path  is  defined  as 

Af(/3)  =  {(x,y,z)  :  {x,y,z)  €  T,  \\Xz  -  ye\\2  <  fiy},  (4) 

where  /3  G  (0, 1)  and  y  =  xT z/n. 

A  very  interesting  result  is  Theorem  9  of  McLinden  [13].  In  the  case  of  linear  program¬ 
ming,  it  states  that  the  central  path  intersects  the  solution  set  at  the  analytic-center,  i.e., 
the  central-path  point  (x(y),  y(y),  z(y))  converges  to  the  analytic  center  (x*,  y*,  z*)  as  y  con¬ 
verges  to  zero.  See  also  Proposition  8.2  in  Megiddo  [15]  and  the  discussion  preceding  it.  This 
fact  plays  a  critical  role  in  the  development  of  most  primal-dual  interior-point  algorithms 
which  attempt  to  follow  the  central  path. 
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3  Algorithm 


The  interior-point  algorithm  proposed  in  this  work  has  been  designed  for  effectively  comput¬ 
ing  the  analytic-center  of  the  solution  set  in  linear  programming  and  it  is  a  modification  of 
the  Kojima-Mizuno-Yoshise  algorithm. 

The  following  lemma  can  be  found  in  Zhang  and  Tapia  [24].  It  provides  a  sufficient 
condition  for  a  strictly  feasible  sequence  {(xfc,  yk,  zk)}  to  converge  to  the  analytic-center  of 
the  solution  set. 


Lemma  3.1  (Zhang- Tapia)  Let  {(xk,yk,zk)}  be  a  sequence  of  strictly  feasible  points.  Let 

(al) 

yk  =  (xk)Tzk/n  -»  0  (5) 

and 

(a2) 

II Xtzk/i,t  -  e||2  0,  (6) 

with  e  =  (1, . . . ,  1)T  e  Rn. 

Then  {(xk  ,yk , zk)}  converges  to  the  analytic-center  of  the  solution  set. 

Our  approach  to  constructing  an  algorithm  for  computing  the  analytic  center  is  to  at¬ 
tempt  to  enforce  conditions  (5)  and  (6). 

A  key  ingredient  in  the  proposed  algorithm  is  the  use  of  the  damped  Newton  method 
applied  to  the  perturbed  Karush-Kuhn- Tucker  conditions 

A  x  -  I,  N 

ATy  +  z  —  c  =0, 

XZe  —  ye  ) 

with  y  >  0.  A  line-search  technique  is  used  as  a  globalization  of  Newton’s  method.  For  fixed 
y  >  0,  the  merit  function  considered  is  //t  :  R2n+m  — >  R  defined  by 

fv{x,y,z)  =  \\Fv{x,y,z)/y\\l- 
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Observe  that 


(7) 


f»(x,y,z)  =  \\{Xz  -  fie)/fi\\l,  for  (x,y,z)  €  T. 

The  function  /M  measures  proximity  or  closeness  to  the  central  path  at  ft. 

The  basic  idea  of  the  algorithm  is  to  consider  a  sequence  of  gradually  shrinking  neighbor¬ 
hoods  {Af(0j)}  of  the  central  path  with  0j  — >  0.  In  the  next  section  we  will  show  that  under 
certain  assumptions  the  proposed  algorithm  generates  an  iteration  sequence  {(xk1yk1zk)} 
such  that  a  subsequence  (xk\  yk] ,  zk>)  belongs  to  M(03).  This  subsequence  is  obtained  by 
considering  a  sequence  of  p's  and  for  each  fixed  p  using  a  damped  Newton’s  method  (with 
a  line-search  globalization  strategy)  to  approximately  solve  the  system  Fn(x,y,z)  =  0. 

A  description  of  the  algorithm  is  the  following: 

Algorithm  3.1  (Long-Step  Shrinking-Neighborhood  (LSSN)  Algorithm)  Consider 
a  strictly  feasible  point  w°  =  (x°,  y°,  z°),  a0  €  (0, 1),  and  0°  €  (0, 1). 

Choose  if  €  (0, 1/2)  and  0  <  /  <  u  <  1. 

Do  until  convergence 

(0)  Set  k  =  0. 

(1)  Define  parameters:  Set  ft  =  pk  =  a°(xk)T zk /n  and  0  =  0k. 

(2)  Check  proximity  to  the  central  path:  If 

\\XkZke/ft-e\\2<0  (8) 


go  to  (6). 

(3)  Compute  new  iterate: 


(3.1)  Solve  the  following  system  for  Xwk  ■ 

=  ( Ark,A/,Azk ): 

*  Ax  ^ 

(  °\ 

F'(xk,  yk,  zk) 

Ay 

=  ~F(xk,  yk,  zk)  +  ft 

0 

V  J 

V e  / 
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(3.2)  Choose  rk  £  (0, 1)  and  compute  the  step-length  ak  =  min(l,  Tkak),  where 

~k  _ zl _ 

min((Xfe)_1  Axk,  (Zk)~1Azk) 

(3.3)  Form  the  iterate 

wk+1  =  (xk+\yk+\zk+1)  =  (xk,ijk,zk)  +  ak{Axk,Ayk,Azk). 

(4)  Line  search: 

(H)  U 

M™k+1)  <  Uwk)  +  yakVUwk)T(Axk,Ajk,Azk)T  (9) 

go  to  (5). 

(4-2)  If  not,  reduce  ak  :=  pak ,  with  p  £  [Z,u]  and  form  the  iterate 

wk+1  =  (x  k+1,yk+1,zk+1)  =  (xk,yk,zk)-hak(Axk,A/,Azk). 

(4-3)  Go  to  step  (4-1). 

(5)  Set  k  :=  k  +  1  ,  flk  =  /3,  pk  —  //.  and  go  to  (2). 

(6)  Set  p  =  fik  =  a°(xk)Tzk /n. 

(7)  Do  (3.1),  (3.2),  (3.3). 

(8)  Decrease  (3 -neighborhood:  Choose  (3k+l  <  (3k . 

(9)  Set  k  :=  k  +  1  and  go  to  (1). 

We  note  that  while  condition  (8)  is  not  satisfied  substeps  (3)  and  (4)  of  the  above  algo¬ 
rithm  are  performed  iteratively  with  the  same  value  of  //,.  It  is  convenient  to  denote  by  I(//,) 
the  corresponding  set  of  indices  i.e., 

IM  =  {£  :  /  =  /<}  •  (10) 

In  the  next  section  we  will  show  that  the  index  set  Z(/j)  is  finite  for  each  //  chosen  by 
Algorithm  3.1.  This  means  that  under  certain  assumptions  the  iteration  sequence  generated 
by  the  LSSN  algorithm,  when  p  >  0  is  fixed,  approaches  the  central  path  in  the  sense  that 
condition  (8)  (or  proximity  to  the  central  path)  is  satisfied. 
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4  Convergence  Properties 

In  this  section  we  prove  that  under  certain  assumptions  the  LSSN  algorithm  generates  an 
iteration  sequence  which  contains  a  subsequence  that  converges  to  the  analytic-center  of  the 
solution  set.  The  convergence  result  is  established  by  showing  that  under  certain  assumptions 
there  exists  a  subsequence  satisfying  conditions  (5)  and  (6)  of  Lemma  3.1.  Therefore  we 
need  to  study  the  behavior  of  the  proposed  algorithm  with  respect  to  the  two  objectives  of 
decreasing  the  gap  and  centering  the  iterates.  In  the  following  two  propositions  we  study 
these  issues.  In  the  first  proposition  we  prove  that  the  gap  sequence  generated  by  the  LSSN 
algorithm  is  a  decreasing  sequence.  In  the  second  proposition  we  prove  that  for  each  fixed 
H  >  0,  the  search  direction  is  a  descent  direction  for  the  function  flt  and  the  sequence 
{fp(xk,yk,zk)}  generated  by  the  algorithm  is  a  strictly  decreasing  sequence.  Also,  we  show 
that  f^{xk,yk,zk)  is  an  upper  bound  for  the  Euclidean  measure  ]|  ~  e|||,  therefore 

condition  (6)  may  be  established  for  any  sequence  of  iterates  {(xfc,  yk ,  zk )}  with  the  property 
that  the  corresponding  sequence  {fP(xk,  yk,  zk)}  converges  to  zero. 

Proposition  4.1  Let  {(xk,yk,zk)}  and  {(/}  be  generated  by  the  LSSN  Algorithm  (Algo¬ 
rithm  3.1).  Then, 

(xk+l)T zk+l  =  (1  —  ok(l  —  ak))(xk)T zk ,  for  all  k>  0  and  for  some  ak  €  (0, 1]. 

Proof 

Define  the  following  sets 

Icp  -  {k  >  0  :  ( xk,yk,zk )  satisfies  condition  (8)  in  the  algorithm  },i.e. 

Icp  =  {k>  0  :  || XkZke/,i  -  e\\2  <  fik),  and  . 

1+  =  {k  >  1  :  (xk~1,yk~1,zk~1)  satisfies  condition  (8)  in  the  algorithm  }  U  {0}. 

For  each  k  we  denote  by  kfp  the  largest  element  in  /+  such  that  kfp  <  k  and  kcp  the 
smallest  element  in  Icp  such  that  k  <  kcp  . 

Let  ak  for  all  k  >  0  be  defined  as 

<7* = a°/  n  (i  -  «m(i  -  for  k?ICP  U  CP 


10 


and 


ak  =  a°,  for  k  e  ICp  U  /+. 

A  fairly  lengthy  but  straightforward  inductive  proof  can  be  used  to  establish  the  following 
results  (a  detailed  proof  can  be  found  in  Gonzalez- Lima  [9]) 

(i)  crk  €  (0, 1],  for  all  k  >  0. 

(ii)  ak  =  1  if  ctk~l  =  1,  for  k  £  Icp  U  /+. 

(iii)  crk+p+1  <  akcp+2  <  ...  <  <  1. 

From  step  (3.1)  in  the  algorithm  we  have 

(. xk+1)Tzk+1  =  (xk)Tzk(  1  -  ofc)  +  no*/i,  for  all  Jfc  >  0.  (11) 

Then,  the  remainder  of  the  proof  consists  in  showing  that 

H  =  ak(xk)Tzk/n,  for  all  k  >  0.  (12) 

If  k  €  Icp  U  /+  then  ak  —  cr°,  and  by  the  construction  of  the  algorithm  we  have 

0(xk)Tzk  k(xkfzk 

a  =  a  -  =  (T  - . 

n  n 

If  k  ^  Icp  U  /^,  assume  without  loss  of  generality  that  k+p  =  0  and  k  =  2. 

From  (11)  we  get 

(x2)Tz2  =  (x1)T21(1  —  a1)  +  (13) 


Multiplying  both  sides  of  equation  (13)  by  (J2/n,  using  the  definition  of  cr2  and  the  fact 


that 


we  obtain 


p.  =  a 


fj,  =  a‘ 


Ax')Tz\ 

n 

(x2)Tz2 

n 


The  use  of  the  same  argument  in  an  inductive  manner  for  k  <  kcp  completes  the  proof  of 

(12).  D 
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Proof. 

Since  wk  €  we  have 

/,(«<*) = -4ii^('”‘)in 

r 

fl 

Therefore, 

v/f(»i)w = -  = -4rc-<“‘)ie  <  ° 

/x^  /xz 

and  (i)  follows. 

Part  (ii)  can  be  demonstrated  in  a  straightforward  manner  from  condition  (9)  of  the 
algorithm. 

Now,  let  us  prove  part  (iii).  Because  ^  |  z  e  is  the  orthogonal  projection  of  Xk Zke  over 
the  subspace  generated  by  e  we  have 

XkZke  -  e  <  |A^Zfce  -  fte  =  /*/M(wfc)1/2, 

n  2  1  2 

for  all  fc. 

Hence, 
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By  Proposition  4.1,  //  =  ak(xk)T zk /n,  with  crk  g  (0,1].  Therefore 


XkZke 
( xk)Tzk/n 


Moreover,  if  =  1  the  previous  inequality  is  satisfied  as  an  equality.  □ 

The  next  lemma  is  needed  for  proving  that  under  mild  assumptions,  condition  (8)  from 
the  algorithm  (proximity  to  the  central  path)  is  always  satisfied. 


Lemma  4.2  Let  p,  >  0  be  a  parameter  selected  by  Algorithm  3.1.  and  assume  that 
min (Xk Zke)/xkT zk  >  7/n  for  all  k  g  T(p)  and  some  7  g  (0, 1), 
where  T(p)  is  the  index  set  defined  by  (10).  Then, 

||(Arfc,  Xyk,  A2fc)||2  <  M  <  +00 
for  all  k  g  J(/<)  and  some  M  >  0. 


Proof. 

Suppose  ( xk)Tzk  >  t  for  some  e  >  0  and  for  all  k  >  0.  By  hypothesis,  xkzk  >  7^  ^  — . , 

hence 

xkzk  >  —  for  all  k  >  0  and  i  =  1,  •  •  • ,  n.  (14) 

n 

Let  T0  —  |(x,  y,z)  g  T  :  xTz  <  (x°)rz0}.  This  set  is  bounded  (see  the  proof  of  Lemma  2.1 
in  Zhang  and  Tapia  [24]).  Moreover,  by  Proposition  4.1,  the  iterates  ( xk ,  yk ,  zk )  g  T0 ■  There¬ 
fore,  equation  (14)  implies  that  {xfe}  and  {zk}  are  bounded  away  from  zero  for  all  k  >  0. 

Let  us  denote  wk  =  (xk,iyk,zk)  and  Xwk  =  (A xk,/Xjk,Azk).  It  follows  from  the  bounded¬ 
ness  of  the  sequences  {xfc}  and  {zk}  that  there  exist  constants  C  and  M\  such  that 


and 


ii(*V)ni2 


/ 

'  A 

0 

0 

\ 

-1 

0 

AT 

I 

V 

zk 

0 

xk 

) 

||An^||2<||(FV)r1||2||^)!|2<M. 


(15) 
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Now,  suppose  that  ( x  k)T zk  — >  0.  By  Theorem  3.2  from  Tapia,  Zhang,  and  Ye  [18]  we  have 
that 

||(A tk,  Azk)\\2  <  Cxok  +  C2{xk)rzk  (16) 

for  some  constants  C\ ,  C2  >  0,  and  ak  defined  as  in  Proposition  4.1  (ak  <  1). 

Equation  (16)  implies  that  there  exists  Ni  >  0  such  that 

||(Arfc,Azfc)||2  <  Ni  for  all  k  >  0.  (17) 

Since  we  are  assuming  that  rank(A)  =  m,  the  matrix  AA1  is  invertible  and 

A?/  =  -(AAT)~1AAzk.  (18) 

Therefore,  (17)  and  (18)  imply  that  there  exists  N2  >  0  such  that 

||  Vlh  <  N2  for  all  lb  >  0.  .  (19) 

The  final  result  follows  from  (15),  (17),  and  (19).  □ 

The  next  lemma  states  that  under  some  natural  assumptions  the  index  sets  J (/i)  are 
finite  for  any  /r  >  0  generated  by  the  LSSN  algorithm. 

Lemma  4.3  Suppose  in  the  LSSN  Algorithm  (Algorithm  3.1)  the  parameter  choice  {t^}  has 
been  made  so  that  the  following  condition  is  satisfied: 

(Cl)  T*  >  r  for  all  k  and  some  r  €  (0, 1). 

Now,  assume  that  the  points  generated  by  the  LSSN  Algorithm  satisfy 

(al)  min {XkZke)jxkT zk  >  7 /n  for  all  k  and  some  7  E  (0, 1). 

Then,  the  index  sets  T{p)  defined  by  (10)  are  finite  for  any  value  of  fi  selected  by  the 
algorithm. 

Proof.  Assume  by  contradiction  that  Algorithm  3.1  selects  a  value  of  g  such  that  l(p)  is 
infinite  i.e., 

r(/()  =  {N,N  +  1,JV  +  2, 
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for  some  N  >  0.  Under  this  assumption  we  will  prove  that 


->  0. 


(20) 


where  /M  is  defined  by  (7)  and  re*'  =  (xk,  yk,  zk).  It  is  easily  seen  that  (20)  implies  that 
condition  (8)  will  be  satisfied  for  k  large  enough  which  contradicts  the  assumption  that  I(p) 
is  infinite. 

Let  us  call  ak  =  min(l, rkak),  where 


A  k 

a  = 


min((Xfc)_1  Axk,  [Zk)~l  Nzk)  ’ 


and  rk  G  (0, 1). 


It  is  known  (see  the  proof  of  Theorem  3.1  in  Zhang,  Tapia,  and  Dennis  [25])  that  Condition 
(cl)  and  Assumption  (al)  imply  that  ak  is  bounded  away  from  zero. 

If  ak  =  ak  in  infinitely  many  iterations,  then  by  using  (ii)  of  Proposition  4.2  we  deduce 
immediately  that  (20)  holds.  Now,  suppose  that  ak  <  ak  for  k  sufficiently  large.  It  is  known 
(e.g.  see  Theorem  6.3.3  in  Dennis  and  Schnabel  [4]  and  the  consequent  discussion)  that  the 
backtracking  line  search  used  in  the  LSSN  algorithm  (Algorithm  3.1)  guarantees  that 

Vffl(wk)TAu,k  Q 
2 

where  Awk  =  (Arfc,  lStjk ,  Azk).  According  to  the  proof  of  Proposition  4‘.2  we  can  write 


V/„(»‘)r  =  -VM)  (22) 

and  by  virtue  of  Lemma  4.2,  we  have  that  ||A»A||2  <  M  for  all  k.  Therefore,  equations  (21) 
and  (22)  imply  that  (20)  is  satisfied  in  this  case  as  well.  This  completes  the  proof  of  our 
lemma.  □ 

Using  Lemma  4.3  we  will  show  that  if  the  parameters  /P  that  define  the  neighborhoods 
of  the  central  path  are  chosen  so  that  they  converge  to  zero,  then  the  iteration  sequence 
generated  by  Algorithm  3.1  possesses  a  subsequence  converging  to  the  analytic-center  of  the 
solution  set.  The  next  theorem  formally  states  this  result. 


Theorem  4.1  Let  {(xk,yk,zk)}  be  generated  by  the  LSSN  algorithm  (Algorithm  3.1)  with 
parameter  choices  a0  G  (0,1),  {flk},  and  {rfc}.  Suppose  the  parameter  choices  have  been 
made  so  that  the  follovhng  two  conditions  are  satisfied: 
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(cl)  Tk  >  r  >  0  for  all  k  and  some  r  £  (0, 1). 

(c2)  0. 

Assume 

(al)  min (XkZke)/xkT zk  >7 /n  /or  all  k  and  some  7  £  (0, 1). 

Then,  there  exists  a  subsequence  of  {(xk ,yk ,  zk)}  that  converges  to  the  analytic- center  of  the 
solution  set. 


Proof. 

By  Assumption  (al),  Condition  (c.1),  and  Lemma  4.3,  there  exists  a  subsequence  {k3} 
such  that  (8)  is  satisfied  at  iteration  k3  i.e., 

,.Xk>Zk>e 


yk> 


-e||2  </?*'. 


According  to  (iii)  of  Proposition  (4.2)  we  have, 

H  Xk’Zk'c 

(xki)Tzki  /n 

and  by  virtue  of  Condition  (c2)  it  follows  that 

Xk*Zk>e 


eh  <  (/e(^J))1/2, 


el|2  — *■  0. 


(23) 


(xki)T  zki  /n 

Now,  for  each  j  >  0,  Proposition  (4.1)  implies  that  (xfcj+1  )T zkj+1  <  (xkJ+1)T zkj+1  and 
(xkj+1)T zkj+1  <  ((xkj )Tzk]  for  some  (  £  (0,1).  Therefore, 

(a :ki)Tzk>  0.  (24) 

The  remainder  of  the  proof  follows  from  (23),  (24),  and  Lemma  3.1.  □ 


5  Complexity  Bound 

In  this  section  we  address  two  important  issues  concerning  the  LSSN  algorithm;  namely  the 
centrality  assumption  in  Theorem  4.1  and  the  complexity  of  the  algorithm.  We  prove  that 
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a  modification  of  the  line-search  strategy  for  the  LSSN  algorithm  leads  to  a  modified  LSSN 
algorithm  which  possesses  polynomial  complexity  and  for  which  the  centrality  assumption 
of  Theorem  4.1  is  always  satisfied. 

In  order  to  facilitate  the  exposition  below  let 


—  {( x,y,z )  strictly  feasible 


min  (XZe) 
xT zjn 


>  7},  for  7  e  (0,1). 


Clearly  if  (xk,yk,zk)  6  IF-,,  for  all  k  >  0,  then  the  sequence  {(xk,  yk,  zk)}  satisfies  Assump¬ 
tion  (al)  (the  centrality  assumption)  of  Theorem  4.1. 

For  each  y  >  0,  the  modified  LSSN  algorithm  presented  in  this  section  considers  an  exact 
line-search  instead  of  the  backtracking  strategy  performed  by  the  LSSN  algorithm.  The  exact 
line-search  is  defined  so  that  the  iterates  generated  by  the  modified  LSSN  algorithm  always 
belong  to  the  set  Therefore,  we  demonstrate  that  the  assumption  of  Theorem  4.1  can 
always  be  enforced.  Moreover,  since  the  merit  function  for  the  line-search  is  the  Euclidean 
norm  proximity  measure  given  by  /#l,  a.  polynomial  bound  can  be  derived  for  the  number  of 
iterations  required  by  the  algorithm  to  satisfy  the  proximity  condition  to  the  central  path 
given  by  (8).  Hence,  polynomial  complexity  is  attainable. 

The  following  is  a  description  of  the  modified  LSSN  algorithm. 


Algorithm  5.2  (Modified  Long-Step  Shrinking-Neighborhood  (LSSN)  Algorithm) 

Consider  a  strictly  feasible  point,  w°  =  (x°,y°,z0),  cr°  £  (0,1)  and  [3°  €  (0,1).  Choose 
7  €  (0,  2  -  y/2]  such  that  7  <  • 

Do  until  convergence 

(0)  Set  k  =  0. 

(1)  Set  yk  =  a°(xk)Tzk /n. 

(2)  Set  y  =  yk  and  (3  —  fdk . 

(3)  Denote  wk’°  =  wk . 

(4)  Set  1  =  0. 

(5)  Set  w  =  (x:y,z)  =  wk’1 . 
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(6)  Check  proximity  to  the  central  path:  If 


go  to  (10). 


\\XZe/p  —  e||2  <  /? 


(7)  Compute  new  iterate: 


(7.1)  Solve  the  following  system  fo 

r  Aiv  — 

*  Ax  ) 

■■  (Ax,  Ay,  Az): 

1  1 

(  0  N 

F\x,y,z) 

Ay 

\  Az  ) 

a. 

+ 

TT 

a 

H 

1 

II 

0 

V  e  / 

(7.2)  Line  search:  Find  the  solution  0  of 


min 


m, 


where  f(B)  =  \\X(9)z(9)/p  -  e||2. 
(7.3)  Form  the  iterate  w+  =  w  +  9 Aw. 

(8)  Set  wk',+l  =  w+  and  6k'‘  =  9. 


(9)  Set  1  =  1  +  1  and  go  to  (5). 

(10)  Set  lk  —  l  and  ivk+1  =  wk,{k . 

(11)  Decrease  0- neighborhood :  Choose  0k  1 1  <  0k. 


(25) 


(26) 


(12)  Set  k  :=  k  +  1  and  go  to  (1). 

Observe  that  in  this  modification  the  iterations  where  the  /3-neighborhoods  of  the  central 
path  are  reached  are  counted  separately  from  the  other  iterations. 

In  order  to  obtain  a  polynomial  complexity  bound  for  the  Modified  LSSN  algorithm  we 
need  to  derive  an  upper  bound  for  the  number  of  iterations  required  by  the  algorithm  to 
satisfy  the  central-path  proximity  condition  given  in  step  (6).  The  idea  for  bounding  the 
number  of  iterations  required  by  the  algorithm  to  approach  the  central  path  is  to  study  the 
line-search  problem  (26)  and  to  find  an  upper  bound  for  its  objective  function  value.  The 
next  lemmas  address  these  issues. 
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Lemma  5.4  Let  7  £  (0,1),  w  —  ( x,y,z )  £  and  g  <  xTz/n. 
Also  let  Aw  —  ( Ax ,  Aj/,  Az)  denote  the  solution  of  the  system 


(  Ax  ) 

( o\ 

F'(x,y,z ) 

Ay 

=  ~F{x,y,z)  +  g 

0 

V  J 

l  e  / 

Consider  <j)  =  ^||A,YAZe||2  where  AX  =  diag(Ax),  A Z  =  diag(Az). 

Then,  w  +  0Aie  £  for  0  <  6  <  0*  =  min(l,  (1  —  7)/^). 

Proof 

Let  in(0)  =  (x(0),  y(0),  z(0))  =  u>  +  0Aio.  Then,  we  have  that 

X(0)z(0)  =  (1  -  0)Xz  +  9ge  +  02AXA z.  (27) 

Since  iv  is  strictly  feasible  AxTAz  =  0  and  from  (27)  we  obtain 

x(0)T  z(0)/n  =  (1  —  0)(x  Tz/n)  +  9g.  (28) 

For  any  0  <  0  <  1  we  deduce  that 

X(9)z{9)  >  7(1  —  0)(xTz/n)e  +  Oge  —  92fge. 

which  implies 

X(0)z(0)  >  7 (x(0)Tz(0)/n)e  +  0[(1  —  7)  —  0<f>]ge.  (29) 

It  follows  from  (29)  that  if  0  <  0  <  0*  =  min(l,  (1  —  7)/^)  then  >  7. 

Now,  let  us  show  that 

x(9)  >  0 ,  z(0)  >  0,  for  all  0  <  0  <  0*  =  min(l,(l  —  7)/^).  (30) 

For  0  =  0  we  have  u>(0)  =  w  and  (30)  is  satisfied  because  of  the  hypothesis  that  iv  £  •TV 
Let  us  suppose  by  contradiction  that  (30)  is  not  satisfied  for  all  9  £  [0,0*].  Then  there 
exists  9  £  [0, 0*]  and  i  £  {1, 2, . . . ,  n}  such  that  either  x;(0)  =  0  or  z;(0)  =  0.  But  then  from 
(29)  we  get 

0  >  7 (x(0)Tz(O)/n)  +  0[(1  -  7)  -  9<ft\g.  (31) 
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From  (28)  we  deduce  that  x(9)Tz(9)/n  >  9 p  so  that  (31)  gives  the  contradiction  0  >  7 Op. 

□ 

Lemma  5.4  states  that  the  minimization  problem  (20)  gives  a  lower  value  of  the  objective 
function  than  does  the  problem 


min  f(9). 


(32) 


The  next  lemma  uses  the  result  obtained  from  Lemma  5.4  and  the  optimization  problem 
(32)  in  order  to  obtain  an  upper  bound  for  the  optimal  value  of  (26).  This  upper  bound  will 
be  used  later  for  bounding  the  number  of  iterations  needed  to  satisfy  (25),  i.e.  the  number 
of  iterations  needed  for  entering  the  /3-neighborhood  of  the  central  path. 

Lemma  5.5  Let  7  £  (0,2  —  \/2] ,  w  =  (x,y,z)  £  p  <  xT z/n  and  p  =  || (Xz  —  pe)/p\\2. 
The  optimal  value  of  problem  (26)  is  bounded  above  by 

^7(1  -7)1 


P  =  pmaxil  -  - 


7 


2x72(1-7)’ 


1  - 


Proof 

From  (27)  we  deduce  that 


Xz  0AXAz 

m  =  \\(l  -9)—  +  9e  +  92 - 

p  p 

11/  /lWXz  oAXAz.. 

=  ll(l-9)(----«)  +  9J— —  lb- 

/X  /X 


-e  2 


As  in  the  previous  lemma,  let  us  denote  <f>  —  ^||AXAZe||2.  Then  for  any  0  <  9  <  1  we 
obtain  the  inequality 

(33) 


f{9)  <  (1  -  9)\\—  -  e\\2  +  92<)>. 

9 

Let  p  —  || (Xz  —  pe)/p\\2  and  p(9)  =  (1  —  9)p  +  92<j>.  From  (33)  we  have 

m  <  p(9). 

Let 

p+  =  min  p(9). 

0  <0<6* 


(34) 


(35) 


It  is  clear  from  (34)  and  Lemma  5.4  that  the  the  optimal  value  of  problem  (26)  is  bounded 
above  by  p+ .  The  remainder  of  the  proof  consists  in  finding  an  upper  bound  to  p+ . 
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Note  that  />(0)  =  /( 0)  =  p.  The  unconstrained  minimum  of  p{9)  =  (1  —  9)p  +  92<f)  is 
obtained  for  9  —  p/{2<f>).  We  consider  four  cases. 

(i)  If  p  <  2(1  —  7)  and  (1  —7 )/<f>  <  1  then  p/(2<f>)  =  9  <  9*  =  (1  —  7)/^  and  consequently 


?  =  PW  =  I1  “ 


(36) 


We  want  to  find  an  upper  bound  for  p+  in  terms  of  p  and  7.  To  this  effect  we  use  the 
inequality 

IIAXAzih  <  JL^XZr'I’iXz  -  ,.e)||s, 

whose  proof  can  be  found  in  [16].  According  to  the  hypothesis  of  our  lemma  w  =  ( x,y,z )  € 
and  p  <  xT z/n  which  leads  us  to  the  inequality 


<!>< 


2x/27' 


(37) 


Substituting  into  (36)  and  using  the  fact  that  p  <  2(1  —  7)  we  obtain  the  following  bound 
for  case  (i): 


f>+  <  (1  - 


7 


(38) 


2V2(1  -7)' 

(ii)  On  the  other  hand  if  p  <  2(1  —  7)  and  (1  —  7)/^  >  1  then  9*  —  1  and  we  may  have 
either  9  <  9*  =  1  in  which  case  (36)  and  (38)  hold,  or  9  >  9*  =  1  and  then 


p+  =  ^(1)  =  <f>  <  p/2. 


(39) 


Note  that  since  7  <  2  —  \/2  then  p/ 2  is  less  than  the  quantity  on  the  right-hand  side  of 
(38),  so  that  (38)  holds  whenever  p  <  2(1  —  7). 

(iii)  Now  if  we  suppose  that  p  >  2(1  —  7)  and  (1  —  '))/</>  <  1,  then  p/{2(j))  =  9  >  9*  = 
(1  —  7)/^,  which  implies 


p+  =  p(9*)  =  (1  - 
Using  (37)  we  obtain  the  upper  bound 


I  ~  7  ,  (1-7) 

<t> 


+ 


1  -  7  > 


,+  <  a  - 


(40) 


(41) 


21 


(iv)  Finally,  if  we  assume  that  p  >  2(1  —  7)  and  (l  —  j)/</>  >  1,  then  p/(2<f>)  =  6  >  0*  —  1, 
and  therefore  (39)  holds  in  this  case  as  well.  But  then  so  does  (38)  and  (41). 

In  conclusion,  we  have  shown  that  (38)  holds  whenever  p  <  2(1  —  7),  while  (41)  holds 
whenever  p  >  2(1  —  7).  We  end  the  proof  of  our  lemma  by  noting  that  the  right-hand  side 
of  (41)  is  greater  than  the  right-hand  side  of  (38)  if  and  only  if  p  >  2(1  —  7). 

□ 

Lemma  5.6  Let  {(xk'l,yk’l,zk'1)}  and  {pk}  be  generated  by  the  Modified  LSSN  algorithm 
with  the  parameter  choice  {/}. 

For  each  k  >  0  let  lk  denote  the  iterate  such  that 

|| Xk',kzk',k/pk  -  e||2  <  /. 

Then  lk  is  well  defined  with 

l 0  =  0(n2  /~f)  and  lk  —  0(n/ 7)  -f  0(log/)  for  k  >  1. 


Proof 

Let  pk’1  =  || Xk,lzk’1  /p,k  —  e || 2 -  From  Lemma  5.5  we  have  that 


and 


JM+i 


2\/2(l  -  7) 


< 


)/■'  if  />*''<  2(1  -7) 
(‘  -  f’kJ  if  />*■'  >  W  -  7)-- 


Using  (43)  and  the  fact  that  /J+1  <  />i  we  0Ltain  the  relation 


/’,+1  < 

It  follows  that 

/’“  <  P 

Therefore,  /’*  <2(1—7)  i11  at  most 

P  - 


'  k,i  _  vm  ~ 

p  pk,  0 


M  ^  nfc,o  _  ,^7(1-7) 


*_  (/’°-2(l-7))/’° 


/27(1  -  7) 


(42) 


(43) 


(44) 


(45) 
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steps. 

From  (42) 


P*"  <  (1  -  2^J_t)  <  (1  -  «‘)'2(1  ~  7) 


where  a  =  7/(2\/2(l  —  7)). 

It  follows  that  pk'‘k+l  <  /3k  in  at  most 

w> 

log(l  -  a) 

steps. 

The  following  observations  are  in  order.  If  pk,°  <2(1  —  7)  <  then  lk  =  lk  —  lk  =  0.  If 
2(1  —  7)  >  f3k  then  P  =  0  and  P  is  given  by  (46). 

On  the  other  hand,  if  pk,°  >2(1—7)  and  2(1  —  7)  <  f3k  we  have  that  lk  =  P  as  defined 
by  (45)  and  P  =  0. 

Hence,  lk  is  well  defined  and  from  (45)  and  (46)  we  obtain  that  lk  <  P  +  P. 

The  rest  of  the  proof  is  devoted  to  finding  bounds  for  pk,°. 

For  any  0  <  p  <  ( x°)Tz°/n ,  it  is  clear  that 

||XV-HI^<II^0||^<((.t0)V)2.  (47) 

Let  p°  =  a°(x°)T z° /n.  It  follows  from  (47)  that  p 0,0  <  n/a°. 

According  to  (45),  this  implies  that 

P  =  0{n2h).  (48) 

If  /3°  is  independent  of  n,  say  /?°  =  1,  then 

l°  =  0(  1)  (49) 


so  that 


Now, 


/°  =  0(n2h). 


(/’°)2  =  \\Xkzk/pk  -  e||i  =  ^y\\Xkzk  -  pke\\\ 


23 


(||X‘z*  -  (*‘V/»)e ||I  +  (/  -  z“  zk /nf  Ml). 


kT  k  . 


By  definition  of  the  Modified  LSSN  algorithm, 


|| Xkzk  -  fik  1e||2  <  I3k  V  \  for  each  k  >  1, 


(51) 


and  because  ( xkT zk/n)e  is  the  orthogonal  projection  of  Xkzk  over  the  subspace  generated 
by  e  we  have, 

|| Xkzk  -  (xkTzk/n)e ||2  <  || Xkzk  -  /"‘elh-  (52) 

Since  pk  =  o°xkT zk /n  we  obtain 

(/•°)2  <  +  n(/  -  x‘V/n)2) 


It  follows  that 


=  7-Jv((Pk  l)2(xkTzk/n)2  +  n(l  -  a°)2{xkT zk /n)2). 

(/'  ) 


k,0\2  ^  (r)2  +  (1  -  o2 


(/t  < 


n 


(53) 


where  /3 *  =  maXA,>0  @k- 

Finally,  from  (45),  (46),  (53),  and  lk  <  lk  +  lk  we  deduce  that  lk  =  0(n/7)  +  0(log/?fc), 

for  A:  >  1.  n 

The  previous  lemma  gives  an  upper  bound  for  the  number  of  iterations  required  by  the 
Modified  LSSN  algorithm  to  enter  the  ^-neighborhood  of  the  central  path.  The  polynomial 
complexity  bound  for  the  algorithm  is  obtained  by  assuming  a  certain  rate  of  decrease  for 
/3k  and  by  finding  a  bound  on  the  number  of  iterations  required  by  the  algorithm  to  reduce 
the  gap  to  a  given  tolerance.  This  is  done  in  the  next  theorem. 


Theorem  5.2  Let  {(.t k,ijk,zk)}  and  {pk}  be  generated  by  the  Modified,  LSSN  algorithm  with 
{(x°,  y°,  2°)}  a  strictly  feasible  point  and  the  parameter  choices  {fik}- 

Choose  e  >  0  and  er°  £  (0, 1)  such  that  cr°((3*  +  1)  <  1  where  (3*  =  maxfc>0  j3k .  Then, 

(i)  xkT zk  <  t  in  at  most  k(  steps,  with  kt  —  0(log(.r°)Tz°/e). 
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(ii)  Let  f3k  be  chosen  as  (3k  =  ((3°)k .  Then,  the  Modified  LSSN  algorithm  iteration  com¬ 
plexity  is 

(x°)Tz°  (x°)Tz° 

0(n2h)  +  0(\nh  +  log  log 


(in)  Let  (3k  be  chosen  as  fik  =  ^jj+T)  ■  Then,  the  Modified  LSSN  algorithm  iteration  com¬ 
plexity  is 

<x°)Tz°  (x°)T-z° 

0(n2/i)  +  0([n/ 7  +  log(log  - - )]  log(- - ). 


Proof 

From  (51)  we  have  that  xk+1zk+1  <  (1  +  fik)fik  =  (1  +  (3k)a°xkT zk /n.  for  all  k  >  0. 
Therefore,  xkT zk  <  gk(x°)Tz°  where  rj  =  (1  +  /?*)tr0  with  =  ma xk>o /3k  and  g  <  1  by 
hypothesis.  Hence,  xkT zk  <  e  in  at  most  ke  =  0(log(x°)r z° /e). 

From  Lemma  5.6 


k<  r.  /r0\ T  O  (  0\X  0 

g/fc  =  o(-iog(i^i-^^ .  > 

k=i  T  6 


0 \k 


-))  +  0(\ogC-^—Y)  if  Pk  =  ((3°) 


(54) 


and 


ke  „  (r°\Tz0  (r°)Tz°  ( X°)T  Z°  ,  1 

£/*  =  0(-  log(^-))  +  0(log(^li-)  log(log(^))  if  fik  =  (55) 

c  t  c 


k= i 


2(k  +  iy 


Items  (ii)  and  (Hi)  follow  from  (54),  (55)  and  the  fact  that  1°  -  0(n2  / 7)  (see  Lemma  5.5). 

□ 


Corollary  5.1  Consider  problem  (1)  and  assume  the  entries  of  the  data  A,b,c  are  integers 
with  L  being  the  bit  length  of  that  input  data. 

Let  {(x°,  y°,  z0)}  be  a  strictly  feasible  point  such  that  (x°)T z°  <  22L . 

Then,  iteration  complexity  of  the  Modified  LSSN  algorithm,  is 

0(n2/ 7)  +  0([n/7  +  L\L),  if  (3k  =  C (3k~l  for  k  >  1  and  for  some  C  €  (0, 1) 

and 

0(n!/7)  +  0([n/7  +  log(I)li).  i//9‘  =  0(±)/or*>  1. 

Proof  Let  e  =  2~2L .  The  proof  follows  directly  from  Theorem  5.2  since  log  -j-  —  <  4 L.  □ 
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6  Numerical  Experience 


In  this  section  we  discuss  the  numerical  results  obtained  from  applying  the  LSSN  algorithm  to 
a  subset  of  the  NETLIB  test  problems.  We  compare  the  performance  of  the  LSSN  algorithm 
with  the  performance  of  the  Kojima-Mizuno-Yoshise  algorithm  and  the  Mizuno-Todd-Ye 
algorithm. 

Our  preliminary  numerical  experimentation  showed  that  the  LSSN  algorithm  as  presented 
in  Section  3  (Algorithm  3.1)  seems  to  perform  better  than  the  Modified  LSSN  algorithm  in 
terms  of  number  of  iterations.  Therefore,  the  results  in  the  present  section  use  the  LSSN 
algorithm  without  modification  (Algorithm  3.1).  At  this  juncture  we  believe  that  the  line- 
search  modification  introduced  in  Section  5  has  more  of  a  theoretical  yalue  than  a  practical 
value.  However,  more  research  will  be  done  on  this  issue. 

6.1  Implementation  Issues 

The  experiments  were  performed  in  64  bit  arithmetic  on  a  Sun  4/model  670-120  using  codes 
implemented  in  MATLAB. 

6.1.1  Stopping  Criteria 

The  algorithm  is  terminated  when 

(\cTxk-bTyk |  || As k-b\\1  \\ATyk  +  zk  -  4U  \\Xkzk  -  ej|A  8 

maX\  1  +  | bTyk\  ’  1  +  IMIi  ’  1  +  ||y*||i  +  Mi’  {xk)Tzk/n  )  ~ 

or  when  the  number  of  iterations  reaches  200. 

6.1.2  Parameter  Choices. 

The  line-search  strategy  (backtracking)  defined  in  step  (4)  of  the  LSSN  algorithm  was  im¬ 
plemented  using  the  value  rj  =  10~4  and  a  fixed  value  p  =  1/2. 

In  all  problems,  the  parameters  rk  were  chosen  as  rk  —  1  —  min{0.05, 0.05(xk)T zk] . 
Therefore  Tk  =  0.95  far  away  from  the  solution  set  and  is  closer  to  1  when  the  duality  gap 
is  small. 
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The  step-lengths  for  the  Mizuno-Todd-Ye  algorithm  were  obtained  following  Ye,  Giiler, 
Tapia,  and  Zhang  [20]. 

The  parameters  fP  were  chosen,  in  general,  as  f3k+1  =  ( (3k )2.  Some  minor  variations,  based 
on  heuristics,  were  considered  in  choosing  /33  when  the  iterates  were  close  to  the  solution 
set.  In  the  numerical  results  presented  in  this  section  /?  =  (3°  =  0.25  for  the  Mizuno-Todd- 
Ye  and  LSSN  algorithms.  The  choice  of  (3°  —  0.25  was  made  so  that  we  could  compare 
the  Mizuno-Todd-Ye  algorithm  with  the  LSSN  algorithm.  Our  experience  has  shown  that 
numerical  efficiency  for  the  LSSN  algorithm  can  be  substantially  improved  by  choosing  (3 
significantly  larger  than  0.25. 

6.1.3  Starting  iterates 

The  starting  points  for  the  LSSN  and  Kojima-Mizuno-Yoshise  algorithms  are  obtained  fol¬ 
lowing  Lustig,  Marsden,  and  Shanno  [12]  and  are  not  necessarily  feasible. 

The  algorithms  generate  a  sequence  of  iterates  that  approach  feasibility  and  have  the 
property  that  the  duality  gap  goes  to  zero.  For  the  LSSN  algorithm,  the  infeasibility  of  the 
starting  iterates  should  be  considered  a  computational  convenience  since  in  all  the  problems 
tested  when  the  initial  /^-neighborhood  of  the  central  path  is  entered  for  the  very  first  time 
the  iterates  are  already  feasible.  Thus,  equivalent  experimentation  could  have  been  done 
starting  the  algorithm  from  the  strictly  feasible  iterates  obtained  in  the  process. 

A  very  interesting  feature  of  the  LSSN  algorithm  is  that  it  provides  a  way  of  computing 
a  point  in  a  /^-neighborhood  of  the  central  path,  for  any  (3  6  (0,1).  In  this  work  we  used 
this  feature  of  the  LSSN  algorithm  for  computing  the  starting  point  for  the  Mizuno-Todd-Ye 
algorithm. 

6.2  Numerical  Results 

Our  numerical  results  are  summarized  in  the  following  three  tables.  Table  1  studies  the 
effect  of  the  choice  of  the  parameter  er°  in  the  performance  of  the  LSSN  algorithm.  Table  2 
compares  the  performance  of  the  LSSN  and  the  Mizuno-Todd-Ye  algorithms.  The  starting 
point  is  the  same  for  both  algorithms  and  the  iterations  are  counted  from  that  point.  Ta¬ 
ble  3  shows  the  performance  of  the  Kojima-Mizuno-Yoshise  algorithm  for  fixed  values  of  the 
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centering  parameters  crk. 

Observe  that  the  worst  performance  of  the  LSSN  algorithm  is  obtained  for  very  small  or 
very  large  values  of  a0.  On  one  hand,  if  the  parameter  cr°  is  small  (e.g.  a0  <  0.0001  for  most 
of  the  problems  we  considered)  the  number  of  iterations  exceeds  the  allowable  maximum 
number  of  200  iterations.  This  large  number  of  iterations  is  due  to  the  fact  that  the  iterates 
approach  the  solution  set  very  fast  but  are  still  far  away  from  the  central  path.  Because  of 
the  proximity  of  the  iterates  to  the  solution  set  the  Jacobian  matrix  is  very  ill-conditioned 
and  the  step-lengths  become  very  small.  On  the  other  hand,  comparing  the  values  of  er°  for 
which  the  algorithm  converges,  it  can  be  observed  that  the  closer  the  value  of  the  parameter 
(T°  is  to  1,  the  stronger  is  the  adherence  of  the  iterates  to  the  central  path  and  the  slower 
is  the  performance  of  the  algorithm  (as  should  be  expected  from  the  theory).  Observe  from 
Table  1  that  er°  equal  to  0.01  seems  to  be  a  good  choice  to  use  in  implementations. 

Notice  that  the  number  of  iterations  required  to  enter  the  initial  /?— neighborhood  does 
not  seem  to  follow,  in  general,  a  constant  pattern  with  respect  to  the  values  of  cr°.  However, 
when  the  line-search  strategy  is  activated,  the  larger  the  value  of  <r°  the  larger  is  the  number 
of  times  the  step  must  be  decreased  (by  backtracking).  We  believe  this  behavior  reflects  the 
fact  that  the  larger  the  value  of  <t°,  the  more  difficult  it  is  for  the  LSSN  algorithm  to  obtain 
appropiate  decrease  in  the  merit  function,  since  a  higher  priority  is  placed  on  centrality. 

For  most  of  the  problems  tested,  the  Mizuno-Todd-Ye  algorithm  requires  a  larger  number 
of  linear  systems  to  be  solved  to  converge  than  does  the  LSSN  algorithm,  as  can  be  observed 
from  Table  2.  Recall  the  cost  per  iteration  of  the  Mizuno-Todd-Ye  algorithm  is  twice  that 
of  the  LSSN  algorithm.  Actually,  for  many  of  the  problems  the  Mizuno-Todd-Ye  algorithm 
also  requires  a  larger  number  of  iterations  to  converge  than  does  the  LSSN  algorithm.  This 
result  seems  to  be  a  consequence  of  the  long  steps  performed  by  the  LSSN  algorithm.  An 
illustration  of  the  behavior  of  the  Mizuno-Todd-Ye  and  LSSN  algorithms  with  respect  to 
closeness  to  the  central  path  and  the  duality  gap  is  given  by  Figures  1  and  2  for  problem 
SHARE2B. 

Finally,  Table  3  shows  that  the  Kojima-Mizuno-Yoshise  algorithm  seems  not  to  be  con¬ 
verging  for  several  of  the  problems  considered.  This  is  a  significant  disadvantage  when 
compared  with  the  LSSN  and  the  Mizuno-Todd-Ye  algorithms.  The  lack  of  convergence  of 
the  Kojima-Mizuno-Yoshise  algorithm  is  due  to  the  fact  that  the  iterates  do  not  stay  in 
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Table  1:  Study  of  the  effect  of  different  choices  of  the  parameter  cr°  for  the  LSSN  algorithm. 


PROBLEM 


APIRO 

0.00001 

20 

0 

2 

0.0001 

21 

0 

3 

0.001 

22 

0 

3 

0.01 

20 

M 

0 

3 

0.1 

25 

8 

0 

5 

0.5 

65 

8 

0 

18 

BLEND 

0.00001 

NC1 * 3 * 5* 

- 

- 

- 

0.0001 

NC 

- 

- 

- 

0.001 

NO 

- 

- 

- 

0.01 

30 

12 

0 

3 

0.1 

40 

12 

2 

6 

0.3 

0.5 

71 

119 

25 

50 

34 

132 

11 

19 

4 

0 

3 

0 

3 

0 

21 

0 

13 

0 

12 

0 

13 

0 

57 

145  ] 

- 

- 

10 

0 

18 

17 

39 

81 

0.00001 

0.0001 

0.001 

0.01 

0.1 

0.5 


0.00001 

0.0001 

0.001 

0.01 

0.1 

0.5 


0.00001 

0.0001 

0.001 

0.01 

0.1 

0.5 


0.00001 

0.0001 

0.001 

0.01 

0.1 

0.5 


0.00001 

0.0001 

0.001 

0.01 

0.1 

0.5 


0.00001 

0.0001 

0.001 

0.01 

0.1 

0.5 


1  Number  of  Iterations. 

^ Number  of  iterations  taken  to  enter  the  initial  /3-neighhorhood  of  the  central  path. 

3 Number  of  times  the  linesearch  strategy  decreases  the  step  length. 

*  Number  of  times  \x  is  decreased. 

5The  allowable  maximum  number  of  200  iterations  is  reached. 


I 
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Table  2:  Comparison  between  the  Mizuno-Tockl-Ye  and  LSSN  algorithms  in  terms  of  number 
of  linear  systems  solved. 


PROBLEM 

M-T-Y 

LSSN 

- -n 

AFIRO 

16 

11 

0.001 

28 

13 

0.01 

32 

17 

0.1 

BLEND 

52 

18 

0.01 

54 

28 

0.1 

54 

69 

0.5 

SCSD1 

62 

21 

0.01 

70 

73 

0.5 

SHARE2B 

18 

15 

0.00001 

42 

15 

0.001 

58 

21 

0.01 

62 

28 

0.1 

64 

64 

0.5 

SCTAP1 

122 

43 

0.1 

90 

34 

0.01 

LOTFI 

106 

46 

0.1 

SCAGR7 

52 

19 

0.001 

84 

67 

0.5 

SCAGR25 

80 

34 

0.1 

SDSDG 

NC1 

47 

0.1 

*  The  allowable  maximum  number  of  200  iterations  is  reached 


Table  3:  Performance  of  the  Kojima-Mizuno-Yoshise  algorithm. 


PROBLEM 

Number  of  Iterations 

- Z* 

AFIRO 

NC1 

0.3 

39 

0.5 

105 

0.8 

BLEND 

NC 

0.3 

NC 

0.5 

NC 

0.8 

SHARE2B 

NC 

0.3 

NC 

0.5 

NC 

0.8 

LOTFI 

NC 

0.5 

NC 

0.8 

SCAGR7 

39 

0.3 

56 

0.5 

146 

0.8 

SCAGR25 

41 

0.3 

53 

0.5 

138 

0.8 

*The  allowable  maximum  number  of  200  iterations  is  reached. 
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or  return  to  a  /^-neighborhood  of  the  central  path  after  such  a  neighborhood  is  reached. 
The  poor  performance  of  the  Kojima-Mizuno-Yoshise  algorithm  is  a  consequence  of  the  fact 
that  our  stopping  criteria  reflects  convergence  to  the  analytic  center,  and  apparently  this 
algorithm  does  not  in  general  converge  to  the  analytic  center. 

In  the  problems  where  the  Kojima-Mizuno-Yoshise  algorithm  converges,  the  number  of 
iterations  increases  as  the  value  of  the  parameter  ak  increases,  as  should  be  expected.  The 
smallest  number  of  iterations  required  by  the  Kojima-Mizuno-Yoshise  algorithm  to  converge 
is  always  larger  than  the  smallest  number  of  iterations  required  by  the  LSSN  algorithm  to 
converge  for  the  same  problem.  It  is  very  interesting  to  note  that  for  those  problems  where 
the  Kojima-Mizuno-Yoshise  algorithm  converges,  the  line-search  strategy  is  never  activated 
when  applying  the  LSSN  algorithm,  i.e.  they  seemed  to  be  rather  easy  problems.  In  fact, 
we  have  found  out  that  those  problems  share  the  feature  of  having  a  unique  primal  solution. 

7  Conclusions 

In  this  work  a  modification  of  the  Kojima-Mizuno-Yoshise  primal-dual  algorithm  is  proposed. 
The  goals  of  approaching  the  central  path  and  the  solution  set  are  combined  to  design  an 
algorithm  which  effectively  computes  the  analytic-center  solution.  The  iterates  generated  by 
the  Long-Step  Shrinking-Neighborhood  (LSSN)  algorithm  proposed  in  this  research  approach 
the  solution  set  while  approaching  the  central  path.  The  approach  to  the  central  path  is  done 
in  such  a  way  that  we  attempt  to  take  full  Newton  steps  (long  steps)  as  much  as  possible. 
Two  aspects  of  the  algorithm  facilitate  this  objective.  Firstly,  when  the  duality  gap  is  large 
we  choose  /^-neighborhoods  that  are  large.  Secondly,  once  we  reach  a  ^-neighborhood  our 
subsequent  step  is  not  required  to  stay  in  a  /^-neighborhood,  i.e.  is  not  required  to  stay 
near  the  central  path.  However,  as  might  be  expected  from  the  Gonzaga  and  Tapia  [8] 
convergence  theory  for  the  Mizuno-Todd-Ye  algorithm,  near  the  solution  when  the  duality 
gap  is  small  full  Newton  steps  actually  stay  close  to  the  central  path  and  often  fall  in  our 
chosen  /^-neighborhoods;  eventhough  the  /^-neighborhoods  are  shrinking. 

Let  us  now  briefly  make  some  comments  concerning  shrinking  /^-neighborhoods  versus 
fixed  /^-neighborhoods  as  employed  in  the  Mizuno-Todd-ye  algorithm.  It  is  clear  that  far 
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from  the  solution  when  the  duality  gap  is  large  it  is  not  computationally  effective  to  enforce 
excessive  fidelity  to  the  central  path;  hence  we  should  allow  for  large  /^-neighborhoods  in  this 
situation.  Our  convergence  theory  given  in  Section  4  for  the  algorithm  described  in  Section  3 
is  based  on  the  Zhang-Tapia  [24]  sufficient  condition  for  convergence  to  the  analytic  center. 
This  theory  requires  the  satisfaction  of  condition  (a2)  given  in  (6);  which  says  that  for  the 
shrinking  /^-neighborhoods  we  must  have  the  sequence  of  /3's  converging  to  zero.  It  is  also 
worth  noting  that  our  polynomial  bound  complexity  theory  given  in  Section  5  requires  the 
sequence  of  /?’ s  to  converge  to  zero.  Quite  recently  Bonnans  and  Gonzaga  [1]  derived  an 
elegant  theory  giving  conditions  that  are  sufficient  for  convergence  to  the  analytic  center. 
Their  theory  is  less  restrictive  than  the  Zhang-Tapia  [24]  theory  and  does  not  require  the 
sequence  of  /Ts  to  converge  to  zero.  Their  theory  suggests  that  in  our  algorithm,  near  the 
solution  when  the  duality  gap  is  small  /?  can  be  held  fixed.  Very  recently  while  this  paper 
was  in  review  we  performed  some  preliminary  numerical  investigations  with  a  version  of  our 
algorithm  which  left  /3  fixed  near  the  solution.  The  results  were  quite  satisfactory.  Hence 
it  seems  as  if  the  Bonnans-Gonzaga,  theory  could  be  used  to  establish  convergence  to  the 
analytic  center  for  a  version  of  our  algorithm  that  allows  fixed  /?  when  the  duality  gap  is 
small.  Moreover,  the  Bonnans-Potra  [2]  convergence  theory,  which  extends  the  Bonnans- 
Gonzaga  theory  to  infeasible  iterates,  could  be  used  to  analyze  our  algorithm  in  the  case  of 
infeasible  iterates.  Recall,  as  we  previously  mentioned,  that  in  all  our  numerical  experiments 
feasibility  was  established  in  the  first  //-neighborhood  pass  and  the  iterates  remained  feasible 
in  subsequent  //-neighborhood  passes.  Hence,  the  theory  presented  in  Section  4  is  applicable. 

A  very  interesting  feature  of  the  LSSN  algorithm  is  that  it  gives  a  practical  way  of 
finding  a  point  in  a  neighborhood  of  the  central  path.  Hence,  the  LSSN  algorithm  can 
be  used  for  finding  a  starting  point  for  the  Mizuno-Todd-Ye  algorithm  and  comparisons 
of  the  performances  of  the  Generic  Kojima-Mizuno-Yoshise,  Mizuno-Todd-Ye,  and  LSSN 
algorithms  have  been  established. 

Numerical  results  showed  that  the  Generic  Kojima-Mizuno-Yoshise  algorithm  is  not  ad¬ 
equate  for  computing  the  analytic-center  solution  even  if  the  centering  parameters  ak  are 
chosen  close  to  one  and  the  iterates  reach  the  central  path  at  some  iteration.  The  Mizuno- 
Todd-Ye  algorithm  and  the  LSSN  algorithm  are,  in  general,  both  capable  of  computing  the 
analytic-center.  The  LSSN  algorithm  compares  favorably  to  the  Mizuno-Todd-Ye  algorithm 
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regarding  total  cost.  In  this  work  we  proposed  the  use  of  the  LSSN  algorithm  for  finding 
the  starting  point  for  the  Mizuno-Todd-Ye  algorithm,  but  different  approaches  of  obtaining 
this  initial  point  may  lead  to  better  performance  of  the  Mizuno-Todd-Ye  algorithm.  The 
performance  of  the  LSSN  algorithm  strongly  depends  on  the  choice  of  the  parameter  a0. 
The  performance  of  the  Mizuno-Todd-Ye  algorithm  depends  on  the  starting  point  and  does 
not  involve  parametric  choices,  thus,  in  this  sense,  the  Mizuno-Todd-Ye  algorithm  is  less 
complicated  than  the  LSSN  algorithm.  However,  it  was  our  experience  that  while  in  theory 
the  Mizuno-Todd-Ye  algorithm  is  guaranteed  to  converge  to  the  analytic  center,  in  practice 
this  is  not  necessarily  the  case  (obviously  due  to  the  ill-conditioning  encountered  near  the 
solution  set). 
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